AFOSR-TR.  C  5  -  1  05  6 


TEXAS  A&M  UNIVERSITY 
College  Station,  Texas 


Q£>0OO&\QjDO 


NONLINEAR  DYNAMIC  RESPONSE 


COMPOSITE  ROTOR  BLADES 


Annual  Letter 


Prepared  by 


Dr.  John  J.  Engblom 
and 

Dr.  Ozden  0.  Ochoa 


of  the 

Mechanical  Engineering  Department 
Texas  A&M  University 


DT  1C 

F'.ECTEl 
0EC0  9  1985 


Submitted  to  the 

Air  Force  Office  of  Scientific  Research 
United  States  Air  Force 


ME  4786-83-10 


Approved  for  public  release } 
distribution  unli»iied* 


Reproduced  From 
Best  Available  Copy 


Contract  No.  F49620-82-K-0032 
October  1983 

85  12  6  02  7 


'4TE 


■  . 

♦  t 

*> 

* 

g 

• 

,s 

> 

-  .  ' 

i*' 

V‘ 

. 

>; 

1 

* 

NONLINEAR  DYNAMIC'  RESPONSE 

OF 

; 

COMPOSITE  ROTOR  BLADES 

*» 

* 

P 

s 

Annual  Letter 

> 

> 

ft 

X 

. 

ft* 

Prepared  by 

'if 

' 

*» 

1 

Dr.  John  J.  Engblom 

1* 

k 

and 

r 

Dr.  Ozden  0.  Ochoa 

-  . 

1 

of  the 

r, 

1 

.  i  ’ 

Si 

T 

.* 

Mechanical  Engineering  Department 

n / 

"d 

tm  ■ 

J 

Texas  A&M  University 

«k 

c 

1 

9 

'  / 

HV 

, 

Submitted  to  the 

fC 

V 

Air  Force  Office  of  Scientific  Research 

*  '» 

Su 

■N 

United  States  Air  Force 

• 

*v 

• 

ME  4786-83-10  Contract  No.  F49520-82-K-0032 

October  1983 

«•** 

AT?  TO?1’?'?  . . .  • 

« 

lir:  >:  •. '•  ■  ' 

.  1 

i 

r ;  •-  • 

* 

d: -i '  r :  --  : 

K  A  ?  —  f 

**, 

c :n.  - 

v 

f 

<-> 

T= 

i  ■* 

- 

.  •  w  •  . 

YY*'-'-Y-Y-Y 
. -Y-W.'-Y-Y 


HCuRlTV  CL  ATION  THIS  rwn#r  u*r»  tm»no 

REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1  NEPORT  NUMBER  _  12  GCVT  ACCESSION  NO. 

jtFOSR  TR-  ., §>/t/)-/r//cl . 

1  RECIPIENT1*  catalog  number 

_ 

4  TlTLt  tm*4  Subtitle, 

Nonlinear  Dynamic  Response  of  Composite  Rotor 
Blades 

5  tyre  or  REPORT  A  RERIOO  COVERED 

iTTtfSftn-Annual  ICalAsu . 

1  Sept  1982  -  31  Aug  1983 

•  rcreorming  otg.  report  number 

ME  4785-83-10 

f  authors; 

John  J.  Engblom,  Ozden  0.  Ochoa 

»  CONTRACT  OR  GRANT  NUMBERf*; 

F49620-82-K-0032 

»  pertorming  oroanij ation  name  and  aodress 

Texas  ASM  University 

Mechanical  Engineering  Department 

College  Station,  TX  77843 

10.  PROGRAM  ELEMENT  PROJECT.  TASY 

AREA  *  WORl^UNlT  NUMBERS 

-  ji'oy&z 

1!.  CONTROLLING  OrriCE  NAME  AND  ADDRESS  /,/A 

Air  Force  Office  of  Scientific  Research/ A/r 
Building  410  ' 

Bolling  AFB,  DC  20332 

12.  REPORT  DATE 

October  1983 

1*.  NUMBER  Or  PAGES 

37  plus  Appendices 

"U.  monitoring  AGENCY  name  •  AOORESS (II  dif/*ran(  Inn i  Controlllni  Olttcn) 

IS.  SECURITY  CLASS,  (ol  Ihl*  nport) 

Unclassified 

ISA.  DECLASSiriCATION-' DOWNGRADING 
SCHEDULE 

It.  distribution  statement  (oi  u»i»  R«pon) 

Approved  for  Public  Release,  Distribution  Unlimited 

17.  ENSTRlIUTION  STATEMENT  (of  the  •fcalrwl  «N#rfd  in  Block  20,  II  tUUerent  from  Report) 

It.  supplementary  notes 

It.  Rev  BOROS  (Conllnu*  on  nrmr,,  •/<#•  11  n,c,,,ny  mtd  Idmllty  *y  Meek  mMntnr) 

^Composite  Materials  Interlaminar  Sfiear  and  Nbrmal  Stresses’ 

Nonlinear  Dynamic  Response  -’Assumed  ^isplacement^and  Hybrid  Models- 

Damage  Mechanisms  ,  ,  >•  ,-*v 

--finite  Elements  \  *  '''  \J 

Large  displacement  Fbrmulation,  '  y,  M  =  u  -  N 

*6.  ABSTRACT  (Contlnu,  an  rid,  II  n*c*«a*'Y  «r<f  Identity  t,  tied  rnmlnr) 

Summarized  are  research  activities  related  to  Nbnl inear  Dynamic  Response 
of  Composite  Rotor  dtadesC Fundamental  to  the  analysis^ s  the  development  of 
a  continuum  formulation  that  can  accurately  account  for  the  effects  of  inter¬ 
laminar  shear  and  interlaminar  normal  stress  variation  thru-the-thickness  of 
a  laminate.  Technical  highlights  of  the  research  efforts  to  date  are  presented 
for  each  of  the^proposed  tasks J  Namely,' Nonlinear  Displacement  Formulation  for 
Ctxnposite  Media,  Incorporate  Damage  Mechanisms  into  Dynamic  Response  Formulation 
and  Correlation  of  Formulated  Response  Model  with  .Experimental  data. -  Also-*' 

teAZrSf  >-  jC  ‘ 


•v  „  *.  , 


w.--  -  - :-'v . 


V>  V 


20.  included  is  a  list  of  papers  and  abstracts  submitted  for  publication/ 
presentation  as  a  result  of  the  first  year  efforts. 


TABLE  OF  CONTENTS 


Page 

1 


> 


£ 


v 

*r 


i 


■.» 


B 


I.  Overview  and  Summary  1 
II.  Summary  by  Task  2 

II. 1.  Task  I:  Nonlinear  Displacement  Formulation  for  Composite  Media  2 


11. 1.1  Continuum  Formulation 

11. 1.2  Large-Displacement  Formulation 

1 1.1.3  Computer  Implementation 

11. 1.4  Analytical  Verification 


II. 3.  Task  III:  Correlation  of  Formulated  Response  Model  with 
Experimental  Data 


III.  References 
IV.  Tables  and  Figures 
V.  Related  Activities 


\ 

V  vj*  ^0  ) 


2 

5 

7 

10 


1 1.2.  .Task  II:  Incorporate  Damage  Mechanisms  into  Dynamic  Response  11 
Formulation 


13 

14 

16 

36 


VI.  Appendices 


IE*  $!*. 


I.  OVERVIEW  AND  SUMMARY 

Fundamental  to  this  work, is  the  development  of  a  continuum  formulation 
that  can  accurately  account  for  the  effects  of  interlaminar  shear  and 
interlaminar  normal  stress  variation  thru-the-thickness  of  a  laminate. 
Furthermore,  emphasis  is  particularly  on  tapered-twisted  airfoil  geometries 
which  can  be  analytically  represented  as  an  assemblage  of  thin  to  moderately 
thick  finite  elements.  To  achieve  solution  efficiencies,  the  elements 
developed  in  this  work  are  of  the  triangular/quadrilateral  plate  type  as 
opposed  to  solid  type  elements. 

Oh  the  basis  of  these  requirements  and  considering  viable  alternatives, 
three  suitable  continuum  formulations  have  been  developed  and  incorporated 
within  a  finite  element  framework.  These  are  herein  denoted  as  the 
(i)  Higher  Order  Displacement,  (ii)  Modif ied-Kirchhoff  and  (iii)  Hybrid 
Stress  formulations,  respectively.  A  computer  code  ,ias  been  developed  to 
test  the  various  elements  on  t:  basis  of  correlations  with  known  analytical 
and  numerical  solutions.  Line'’.-  s.  'll -displacment  equations  have  been  , 
implemented  in  the  code1 and  many  wests  have  been  performed.  It  is  noted 
that  the  code  has  some  unique  features,  e.g.,  it  can  assemble  elements 
having  an  unequal  number  of  degrees  of  freedom  at  its  nodes,  it  treats 
arbitrary  ply  orientations  and  it  . performs  integration  on  a  layer-by-layer 
basis  through  the  laminate.  Herein  a  layer  refers  to  either  a  lamina  or 
to  a  sub-set  of  laminae  having  equal  p  y  orientations.  The  latter  feature 
is  essential  in  developing  a  fully  nonlinear  capability. 

Significant  efforts  have  also  been  devoted  to  developing  a  suitable 
large  displacement  formulation.  Due  to  the  requirement  that  interlaminar 
stresses  be  accurately  represented,  a  total . Lagrangian  formulation  is 
utilized  and  is  based  upon  the  complete  Green's  strain  tensor.  A  geometric 
and  large-displacement  stiffness  formulation  has  been  implemented  in  f-e 
computer  code  based  upon  a  form  of  the  nonlinear  strain-nodal  displacement 
relationships  suitable  for  each  of  the  elements  under  development. 

An  extensive  literature  sui vey  has  been  performed  to  identify  analytically 
tractable  methods  of  treating  damage  accumulation  in  composites.  Since 
emphasis  in  this  work  is  on  the  development  of  incremental  response  solutions, 
the  computational  approach  must  Have  the  capability  to  (i)  predict  and 
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differentiate  between  relevant  failure  modes,  (ii)  modify  constitutive 
equations  appropriately  and  (iii)  perform  equilibrium  iterations  to  assume 
stress  redistribution  based  upon  the  extent  of  damage.  Use  of  "piecewise 
smooth"  failure  criteria  in  conjunction  with  "damage  state"  variables 
provide  a  good  basis  for  incrementally  tracking  damage.  This  approach  is 
currently  being  formulated  for  incorporation  in  the  computer  code.  Note 
that  integration  for  an  element  is  performed  on  a  layer-by-layer  basis 
which  allows  for  J'mage  effects  to  be  characterized  at  the  layer  level. 

Experimental  data  of  the  type  required  to  substantiate  damage  predictions 
has  been  assembled  to  the  extent  possible.  Analysis/test  correlation 
efforts  will  be  performed  when  the  nonlinear  formulation  including  damage 
effects  is  fully  implemented.  It  is  noted  that  useful  experimental  data  ■ 
is  quite  limited. 

Technical  progress  in  this  program  has  been  substantially  on  schedule. 

It  is  believed  that  the  originally  proposed  three  year  program  can  be 
completed  within  the  given  time  frame. 

II.  SUMMARY  BY  TASK 

This  section  presents  technical  highlights  of  the  research  efforts  to 
date  for  each  of  the  three  tasks.  Details  of  the  analytical  formulation 
are  presented  in  the  Appendices. 

II. 1.  TASK  I:  Nonlinear  Displacement  Formulation  for  Composite  Media 
II. 1.1  Continuum  Formulation 

Two  variational  principles,  the  principle  of  Minimum  Potential  Energy 
and  the  Principle  of  Modified  Complementary  Energy,  are  used  to  develop 
two  distinctly  different  finite  element  models,  the  assumed  displacement 
model  and  the  hybrid  stress  model  respectively.  These  models  incorporate 
the  effects  of  transverse  shear  and  normal  deformations  whose  contributions 
are  recognized  as  essential  for  accurate  laminate  analysis  [1-6]. 

Within  each  formulation,  element  stiffness  and  force  matrices  are 
determined  for  each  element these  matrices  are  then  assembled  to  represent 
the  final  system  cf  equations  and  a  solution  procedure  for  the  unknown 
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nodal  displacements  are  provided.  Coordinate  transformations  to  describe 
ply  orientations  of  a  composite  media  are  taken  into  account.  The  in-plane 
stresses  are  calculated  from  constitutive  relations  of  orthotropic  continuum 
whereas  transverse  shear  and  normal  stresses  are  calculated  from  equilibrium 
considerations. 

The  developed  elements  are  tested  for  linear  static  and  dynamic  analysis. 
The  test  problems  and  the  results  are  presented  in  Section  II. 1.4.  The 
finite  element  models  are  herein  briefly  discussed. 

i.  Assumed  Displacement  Model 

A.  Higher  Order  Displacement  Formulation 

The  thru-the-thickness  effects  can  be  incorported  into  an  analysis  by 
choosing  a  displacement,  field  that  eliminates  two  major  shortcomings  of 
the  classical  plate  theory;  namely  normals  remain  normal  and  in-plane 
displacements  are  linear  thru  the  thickness.  These  shortcomings  are 
eliminated  by  prescribing  independently  the  reference  surface  displacements 
and  rotations  of  the  normal  and  including  higher  order  terms  for  in-plane 
displacements.  This  is  accomplished  by  the  following  variation 

u(x,y,z)  =  uQ(x,y)  +  z.x(x,y)  +  7 ->x(x,y) 
v(x,y,z)  =  v0(x,y)  +  z\,{x,y)  +  z  2iy(x,y) 
w(x ,y,z)  =  Wq ( x ,y ) 

The  neutral  surface  displacements  are  represented  by  u0,  v0  and  w0,  the 
rotation  about  y-axis  is  denoted  by  'i'x  and  the  rotation  about  the  x-axis 
is  *Py .  The  coefficeients  of  z2,  i.e.,  0X  and  ty,  are  contributions  from 
transverse  deformations  [5,6]. 

The  elements  developed  are  designated  as  the  quadrilateral  higher 
order  displacement  (QHD)  models.  QHD40  is  an  eight-noded  element  with 
seven  degrees  of  freedom  (three  midsurface  displacements,  two  rotations 
and  two  higher  order  terms  for  in-plane  displacements)  per  corner  node  and 
three  degrees  of  freedom  (transverse  midsurface  displacement  and  two 
rotations)  per  mid-side  node.  1  Element  QHD28  is  a  simplified  version  of 
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QHD40  where  the  mid-side  nodes  are  eliminated.  It  should  be  noted  that  when 
the  two  higher  order  terms  for  in-plane  displacements  at  each  corner  node  are 
omitted,  QHD28  reduces  to  the  widely  used  four-noded  bilinear  plate  element 
(QHD20) . 

The  transverse  shear  and  normal  stresses  of  QHD40  display  a  cubic 
variation  thru-the-thickness.  The  displacement  field,  nodal  degrees  of 
freedom  and  the  resulting  stress  fields  are  stated  in  Appendix  IA. 

B.  Modified-Kirchhoff  Formulation 

The  Kirchhoff-Love  assumption  for  normals  to  the  reference  surface  is 
relaxed  by  incorporating  shear  rotations  as  additional  degrees  of  freedom 
in  the  formulation  [9].  Thus  the  assumed  displacement  field  allows- the 
transverse  shear  deformations  but  neglects  the  transverse  normal  deformations. 
The  rotations  yx  and  yy  are  incorporated  in  the  displacement  variation  as 
follows . 

w(x,y)  =  w0(x,y) 

u(x,y,z)  =  u0(x,y)  -  z(~  +  >x) 

v(x,y,z)  =  v0(x,y)  -  z(t^  +  yy) 

The  transverse  displacement  w(x,y)  is  chosen  such  that  it  will  guarantee 
plausible  stress  fields  which  will  characterize  the  transverse  effects 
accurately. 

This  approach  is  implemented  in  the  formulation  of  an  eight-node 
quadrilateral .element  with  32  degrees  of  freedom-  QD32,  a  six-node  triangular 
element  with  27  d.o.f.  -  TD27  and  a  seven-node  triangular  element  with 
27  d.o.f..  -  TD27M.  The  stress  fields  obtained  for  these  elements  represents 
a  quadratic  thru  the  thickness  variation  for  the  transverse  shear  stresses  -and  a 
cubic  variation  for  the  transverse  normal  stress.  The  respective  displacement 
fields,  nodal  degrees  of  freedom  and  stress  fields  are  given  in  Appendix  IB. 

ii.  Hybrid  Stress  Model 

In  this  formulation  a  stress  distribution  within  the  interior  of  the 
element  is  expressed  in  terms  of  finite  parameters  such  that  equilibrium 
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,is  satisfied,  also  an  assumed  displacement  distribution  is  used  on  the 
boundary  of  the  element  expressed  in  terms  of  generalized  nodal  displacements 
such  that  the  interelement  compatibility  is  retained  [4]. 

The  element  developed,  QHS32,  is  a  four-node  quadri lateral  with  32 
degrees  of  freedom.  In  addition  to  an  assumed  displacement  field  it  has 
a  26-parameter  stress  field  which  provides  cubic  variation  for  transverse 
shear  stresses  and  a  quartic  variation  for  transverse  normal  stress  through 
the  thickness  of  the  laminate.  The  stress  field  along  with  the  assumed 
displacement  variation  is  stated  in  Appendix  II. 

II. 1.2.  Large  Displacement  Formulation 

Inclusion  of  geometrically  nonlinear  effects  in  the  formulation  must 

be  based  upon  both  the  geometry  to  be  analyzed  and  upon  the  type  of  stress 
prediction  capabilities  desired.  The  classical  approach  to  thin  plate 
analysis  has  been  to  use  the  Kirchhoff-Love  assumptions  in  conjunction 
with  the  nonlinear  von  Karman  relations  [11,12].  As  previously  indicated, 
the  Kirchhoff-Love  assumptions  are  relaxed  in  this  work  to  allow  for  a 
more  accurate  definition  of  interlaminar-shear  and  interlaminar-normal 
stress  variations.  These  stresses  can  vary  substantially  through-the- 
thickness  for  the  geometries,  of  interest,  i.e.,  thin  to  moderately  thick 
plate  type  structures.  Furthermore,  the  requirement  that  these  stresses 
be  accurately  determined  means  that  the  nonlinear  portion  of  the  strain- 
displacement  relationship  must  contain  all  significant  coordinate  displace¬ 
ments.  The  complete  Green's  strain  tensor  is  utilized  in  this  work,  therefore, 
to  account  for  all  significant  contributions  to  the  interlaminar  stress  field. 
With  respect  to  fixed  Cartesian  coordinates  x,  y,  and  z,  the  strain  tensor 
has  the  form 


_  3U  +  3V  h  f  3U  3U  +  3v  3v  ]  3w  3w' 

Yxy  "  ay  ax  lax  ay  ax  ay  ax  ay 

where  u,  v  and  w  represent  displacements  in  the  x,y,z  coordinate  directions, 
respectively.  Mote  that  the  other  strain  components  are  obtained  by  a 
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suitable  permutation.  In  sma il-displacement  analysis,  the  quadratic  terms 
are  neglected  to  give  simply  the  linear  strain  approximation. 

Based  on  the  Green's  strain  tensor,  the  strain  to  nodal  point  displacement 
relationship  can  be  specified  for  elements  under  development.  It  takes  the  form 

{e>  =  [B]{a} 

where  (e}  is  the  vector  of  strain  components,  {A}  the  vector  of  nodal  point 
displacements  and  [B]  a  function  of  derivatives  of  the  element,  shape  functions. 
The  quadratic  terms  in  the  strain  tensor  result  in  [B]  being  a  function  of 
displacement  state  and,  therefore,  an  incremental  equilibrium  formulation  is 
required.  The  incremental  strain-nodal  displacement  relationship  takes  the 
form 

{5c}  =  ([Bq]  +'  [Bl])  {5A} 

where  {us:}  and  { 6 A }  represent  incremental  strains  and  nodal  displacements, 
respectively,  [B  ]  and  [BJ  are  the  small  and  large  displacement  contributions 
to  the  incremental  strains.  Based  on  the  incremental  equilibrium  equations, 
the  displacement  formulation  gives  the  force-displacement  relationships 

[K0]  =  /  [Bq]T  [U][B0]  dV 
v 

[Kl]  =  /  ([B0]T  [DJ[Bl]  +  [BL]T[D][BL]  +  [BL]T[D][Bo])  dV 

V 

where  [D]  is  an  elasticity  matrix  obtained  simply  from  the  constitutive 
equations  and  integration  is  over  the  volume  V  of  the  element.  [Kq]  is 
denoted  the  small-displacement  stiffness  matrix  and  [K^]  is  denoted  the 
large-displacement  stiffness  matrix.  Since  response  is  also  a  function 
of  stress  state,  the  geometrical  stiffness  matrix  [Kg]  is  required  and  is 
obtained  from 

[Kg] { 5A }  = 

v 


/d[BL]T{o}dV 


where  {a} is  the  vector  of  stress  components.  Note  that  the  hybrid  stress 
formu1ation  similarly  gives  force-displacement  forms  involving  the  stress 
and  displacement  state. 

Inertial  effects  are  analytically  trebled  as  a  mass  matrix  [M]  which 
is  a  function  of  density  and  the  element  shape  functions  (see  Appendix  III). 
These  matrix  forms  are  required  in  formulating  static/dynarnic  response 
solutions  and  the  incremental  equilibrium  equations  have  the  general  form 

+  ([kq1  +  [Kl]  +  EKq]) 

where  the  mass  and  stiffness  matrices  represent  an  assembly  of  the  elemental 
matrices  previously  discussed,  { 6 u }  and  { 6 u }  represent  the  incremental 
displacements  and  accelerations  for  the  mathematical  model  and  {SF} 
represents  the  vector  of  incrementally  applied  forces. 

In  developing  a  geometrically  nonlinear  formulation,  the  effort  is 
largely  in  defining  the  incremental  strain-nodal  displacement  relationship. 
Having  developed  this  relationship  for  a  particular  element,  stiffness  matrices 
are  readily  developed  as  the  preceding  equations  indicate.  These  rel utiunshlps 
are  presented  in  Appendix  IV.  The  form  of  these  equations  is  the  same  for 
all  elements. 

II. 1.3.  Computer  Implementation 

A  comruter  code  has  been  developed  for  the  purpose  of  implementing  the 
various  continuum  .formulations .  At  present,  the  code  performs  the  following 
functions: 

(i)  element  stiffness  matrix  generation 

(ii)  element  mass  matrix  generation 

(iii)  assembly  of  equilibrium  equations 

(iv)  decomposition  and  solution  of  equilibrium  equations 

(v)  fundamental  frequency  and  mode  shape  calculation 

A  characteristic  of  the  elements  under  development  is  that  node  points 
can  have  different  numbers  of  degrees  of  freedom,  i.e.,  typically  mid-side 
nodes  have  fewer  degrees  of  freedom  than  corner  nodes.  The  code  has  been 
fashioned  to  handle  this  condition.  All  cf  the  integration  is  performed  on 
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a  1ayer-t>y-layer  basis  thru  the  thickness  of  the  laminate.  This  approach 
is  fundamental  to  developing  the  capability  to  allow  for  inelastic  material 
behavior  and,  ultimately,  to  the  inclusion  of  damage  mechanisms  in  the 
formulation. 

Since  solution  of  the  equilibrium  equations  is  a  vital  component  in 
the  overall  solution  strategy,  it  is  appropriate  to  discuss  the  numerical 
methodology  used  in  solving  these  equations.  The  intent  is  to  obtain  a 
higher  ordered  variation  of  the  transverse  shear  and  normal  stresses 
(axz,  OyZ,  and  o2Z)  than  can  be  obtained  via  the  equilibrium  equations. 

The  solution  procedure  can  be  thought  of  as  described  below.  Assume  that 
the  in-plane  stresses  (oxx,  Cyy,  and  cxy)  within  each  layer  of  a  particular 
element  have  been  determined  at  selected  locations,  i.e.,  through  solution 
of  the  constitutive  equations.  In  the  code  as  presently  written,  these 
locations  are  specified  as  tne  element  centroid  and  element  nodal  points. 
The  equilibrium  equations  (in  the  absence  of  body  forces)  have  the  indicia! 
form 


°ij»j  =  0 

from  which  it  follows  that  the  thru-the-thickness  shear  stress  variation  can 
be  written  in  numerical  form  for  the  it*1  layer  as 

Aoxz-j  =  ~(cxx*x  +  Gxy’y)i  •  ^i 
and 

Acyzi  =  "(axy,x  +  Gyy ^i 

Here,  the  left-hand-side  represents  the  change  ’n  stress  from  the  lower  to 
the  upper  surface  of  the  i^  layer  and  AZ.j  is  the  thickness  of  the  i^  layer 
at  a  particular  location.  The  derivatives  with  respect  to  x  and  y  in  the 
expressions  above  are  readily  computed;  this  is  because  in-plane  stresses 
within  a  layer  are  related  to  element  displacements  through  derivatives  of 
element  shape  functions  in  conjunction  with  a  material  definition. 


9 


For  an  n  layered  laminate,  n  equations  can  be  written  ir  terms  of 
both  the  unknown  shear  stresses  at  layer  interfaces  and  the  shear  stresses 
at  the  laminate  surfaces.  Assuming  the  laminate  has  shear-free  surfaces, 
the  equations  above  give  n  equations  in  n-1  unknowns,  so  that,  the  equation 
set  is  over-determined.  The  equations  have  the  matrix  form  below. 


n  x  (n-1) 


(n-1)  x  1 


(n  x  1) 


where  Ixzi  =-(oxx,x  +  oXy,y).  AZ^  and  ox2 .  represents  the  shear  stress 
acting  at  the  interface  of  the  j-lth  and  layer.  A  similar  equation 
set  is  obtained  by  replacing  oXZj  with  OyZ^  and  Ixz.  with  IyZi-  These 
equations  are  solved  by  utilizing  a  least-squares  orthonormalization 
procedure  [13].  Due  to  the  simplicity  of  the  terms  in  the  coefficient 
matrix,  a  concise  closed-form  solution  is  obtained.  Having  determined  the 
transverse  shear  stresses,  the  transverse  normal  stress  variation  is 
determined  through  the  numerical  form  of  the  third  equilibrium  equation 
for  the  i^h  layer. 


Aazzi  _  ~(°xz’X  +  °yz»y) j 

As  before,'  the  left-hand-side  represents  the  change  in  stress  through  the 
i th  layer.  Appropriate  polynomial  functions  are  utilized  to  describe  the 
oxz  and  OyZ  in-plane  variation.  These  functions  are  differentiated  to 
obtain  the  right-hand-side  of  the  equations  above.  Again  the  equation  set 
is  overdetermined  because  the  normal  tractions  are  known  at  the  laminate 
surfaces.  Solving  for  czz  proceeds,  therefore,  in  identically  the  same 
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manner  as  discussed  in  calculating  rX2  and  0yZ.  Parenthetical ly,  inclusion 
of  body  forces  at  a  later  date  can  be  accomplished  with  little  difficulty. 

II. 1.4.  Analytical  Verification 

An  approach  to  the  successful  application  of  Higher  Order  Displacement 
type  elements,  i.e.,  for  thin  to  moderately  thick  geometries,  is  to  utilize 
reduced  numerical  integration  where  as  this  is  not  necessary  for  Modified- 
Kirchhoff  Formulation.  This  approximation  technique  brings  along  the  choice 
of  implementing  it  overall  or  selectively  to  the  strain  energy  components. 

For  the  QHD  formulation,  only  the  transverse  shear  components  are  integrated 
with  reduced  order.  The  integration  order  may  affect  the  physical  behavior 
of  the  element  by  introducing  spurious  zero  energy  modes.,  It  is  desirable 
to  have  only  rigid  body  modes  since  there  does  not  yet  seem  to  be  a  generally 
accepted  method  of  controlling  the  additional  modes.  For  QHD40,  3x3  Gaussian 
quadrature  along  with  the  2x2  quadrature  for  the  transverse  shear  components 
is  employed.  QHD28  and  QHD20  formulations  are  similarly  integrated  with 
2x2  and  1x1  Gaussian  quadratures.  ,  Manipulation  of  quadrature  rules  may 
result  in  undesirable  element  behavior.  The  presence  of  spurious  zero 
energy  modes  in  addition  to  the  rigid  body  modes  may  detract  from  overall 
performance.  A  spectral(eigenvalue)  test  is  conducted  with  and  without  full 
quadrature  to  observe  the  zero  energy  modes  of  the  QHD  elements.  The 
quadrature  order,  the  number  of  zero  eigenvalues  and  the  corresponding 
number  of  spurious  zero  energy  modes  for  QHD40,  QHD28  and  QHD20  are  listed 
in  Table  1.  The  spurious  modes  of  QHD28  are  illustrated  in  Figure  1. 

It  is  also  noteworthy  to  observe  the  effect  of  reduced  integration  on 
the  representation  of  the  generalized  forces.  In  order  to  illustrate  the 
effect,  the  forces  associated  with  the  transverse  displacement  of  a  corner 
node  are  sketched  in  Figs.  2A  and  2B  for  QHD28  with  and  without  reduced 
integration  respectively. 

To  demonstrate  the  performance  of  QHD40,  QHD28  and  QD32  in  comparison 
to  classical  piate  theory  and  elasticity  solutions,  the  example  problems 
of  Table  2  are  solved. 


>  “ 
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Example  1.  Simply  Supported  Square  Plate  with  Sinusoidal  Load. 

The  results,  given  in  Tables  3  and  4  are  comparisons  of  the  elasticity 
solution  of  Pagano  to  those  calculated  via  the  QHD40,  QD32  formulations. 

Figures  3A-B  and  4A-B  display  the  transverse  normal  stress  variation  for 
QHD40  and  QD32  respectively  at  the  center  and  at  a  point  on  the  edge  of 
the  plate  respectively. 

Example  2.  Simply  Supported  Square  Plate  with  Uniform  Load. 

The  results  presented  in  Figure  5  and  Figure  6  display  the  difference 
between  the  QHD  formulation  and  the  classical  plate  theory.  The  bending 
stresses  displayed  in  Figure  6  for  QHD28  are  in  close  agreement  with  those 
of  QHD40  but  the  transverse  shear  stresses  of  Figure  5  present  rather  poor 
correlations.  The  transverse  shear  stresses  and  bending  stress  of  QHD32 
are  as  shown  in  Figures  7  and  8. 

Example  3.  Cylindrical  Bending 

The  transverse  shear  stress,  aX2,  calculated  via  the  QHD40  and  QD32  . 
formulations  is  illustrated  in  Figure  9  and  10  respectively.  The  bending 
stress  rxx  of  QHD40  is  displayed  in  Figure  11  and  oxx  of  QD32  is  shown  in 
Figure  12.  The  results  are  compared  to  the  elasticity  solution  obtained 
by  Pagano. 

Verification  of  the  dynamics  portion  of  the  computer  code  has  been 
limited  thus  far  due  to  the  emphasis  placed  on  obtaining  good  statical 
results.  A  first  check  has  been  performed,  however,  in  calculating  the 
fundamental  frequency  of  a  simply  supported  (isotropic)  plate.  The 
calculated  frequencies  for  elements  QHD28,  QHD40  and  QD32  are  all  within 
2%  of  the  closed-form  solution  for  a  simple  36  element  plate  model. 

II. 2.  TASK  II:  Incorporate  Damage  Mechanisms  into  Dynamic  Response  Formulation 

The  literature  survey  performed  has  been  quite  helpful  in  terms  of 
delineating  the  viable  approaches  to  including  damage  mechanisms  in  the 
analysis.  Relevant  failure  modes  of  interest  include  those  listed  below. 
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(i)  fiber  fracture 

( i i )  fiber-matrix  debonding 

(iii)  matrix  cracking  (parallel  and  transverse  to  fibers) 

(iv)  delamination 

Several  smooth  failure  criteria,  e.g.,  [14-17]  have  been  developed  in 
recent  years  to  represent  the  failure  of  composites.  These  criteria,  to 
varying  degrees,  can  predict  "failure"  but  do  not  identify  a  particular 
mode  of  failure.  In  performing  incremental  "damage"  analysis,  it  is 
essential  to  both  predict  failure  and  to  characterize  it,  e.g.,  do  fibers 
rupture,  dees  delamination  occur,  etc.  The  computational  approach  must, 
therefore,  differentiate  between  viable  failure  mode1'  and  appropriately 
alter  the  constitutive  equations  on  an  incremental  basis.  This  can  be 
accomplished  by  impementing  a  piecewise  smooth  failure  criteria,  e.g., 

[18]  in  the  finite  element  formulation.  The^ general  failure  criteria  is 
then  comprised  of  m  separate  inequalities  of  the  form 

Fj  ({a})  <1  ;  j  =  1,2,. . .  ,m 

at  each  layer  level  within  each  element.  These  criteria  can  differentiate, 
between  (i)  tensile  and  compressive  fiber  failure,  (ii)  tensile  and  compressive 
matrix  failure  and  (iii)  delamination  at  layer  interfaces. 

As  progressive  damage  occurs  throughout  incremental  loading  (whether 
it  be  static  or  dynamic),  it  is  essential  that  violation  of  failure  criteria 
inequalities  be  reflected  in  modification  of  the  material  properties.  This 
can  be  achieved  by  including  damage  state  variables  [19]  in  the  constitutive 
equations  to  reflect  "stiffness  reduction."  These  equations  can  be' 
represented  as  , 

(o)  =  [D][y]{e} 

where' [D]  represents  the  material  matrix  and  [y]  contains  the  damage  state 
variables.  The  letter  provide  the  basis  for  changing  the  D-jj  terms  based 
upon  the  extent  to  which  the  failure  criteria  are  violated. 
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In,  conjunction  with  the  above  it  is  essential  to  perform  equilibrium 
iterations  within  each  analysis  increment.  This  is  required  to  assure  that 
stress  redistribution  is  properly  accounted  for  as  damage  progresses. 
Currently  the  approach  to  modelling  "damage  mechanisms"  is  under  development 
for  implementation  in  the  formulations. 

II. 3. 3  TASK  III:  Correlation  of  Formulated  Response  Model  with  Experimental 
Data 

Some  quantitative  data  relating  to  the  impact  damage  of  composite 
specimens  has  been  assembled  [21-23],  It  will  be  utilized  along  with  any 
additional  data  obtained, to  perform  analysis/test  correlations.  Since  the 
nonlinear  formulation  including  damage  effects  is  not  complete,  no  use  of 
the  test  data  has  been  made  to  this  point.  Much  effort  has  been  devoted, 
however,  to  correlating  the  analysis  wit'.i  both  closed-form  and  numerical 
solutions  as  previously  discussed. 
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Table  1,  Spurious  Zero  Energy  Nodes  of  the  QHO  Family 


Quadrature  Order 

Number  of  Zero 
Eigenvalues 

Number  of  Spurious 
Modes 

0HD4O 

3x3  with  2x2 
for  transverse 
shear  terms 

6 

0 

QHD28 

2x2  with  lxl 
for  transverse 
shear  terms 

9 

3 

QHD20 


2x2  with  lxl 
for  transverse 
shear  terms 


8 


2 


Quadrant  of  a  plate  .  Quadrant  of  a  plate  Half  a  strip 


ess 


Figure  1.  Spurious  zero  energy  modes  of  QHD28 


Figure  2.  Generalized  forces,  associated  with  the  transverse 

displacement  w0  of  a  corner  node  with  full  quadrature 
.  as  shown  in  (A)  and  with  reduced  integration  as  in 
(B)  for  QHD28.  , 

Figures  3  &  4.  Transverse  normal  stress,  azz  ,  variation  for  QHD40 
and  QD32  respectively,  at  (A):  center  of  the  plate 
(oz  =  oz/100)  and  at  (B):  a  point  on  the  edge  of 
the  plate  (oz  =  10oz). 

Figures  5  &  6.  Thru-the-thickness  oxz  variation  in  element  #6. 

OyZ  variation  in  element  #31;  in  comparison  to 
CPT  results  (oxz  =  ax2/100,  ~rs  =  Oy7/ 100)  for 
0HD40  and  QD32  respectively. 

Figures  7  8  8.  Comparison  of  bending  stresses  of  QHD40,  QHD28, 
and  QD32  to  that  of  CPT  in  element  #36. 

(ox  =  ox/100) 

Figures  9  8  10.  Transverse  shear  stress,  axz,  for  QHD40  8  QD32  in 
comparison  to  elasticity  solution  of  Pagano  in 
.  element  #20.  (cX7  =  oxz/100) 

Figures  11  8  12.  Bending  stress,  oxx,  .variation  of  QHD40  and  QD32 
respectively  in  comparison  to  elasticity  solution 
in  element  #1.  (ax  =  ax/100) 
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The  schedule  of  work  performance  is  presented  below  for  the  three  year 
program.  This  schedule  is  as  originally  proposed  as  no  changes  £re  requested. 
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APPENDIX  IA 


HIGHER  ORDER  DISPLACEMENT  MODELS 


.  QHD40 

NODAL  DEGREES  OF  FREEDOM: 

Corner  Nodes  -  {u0  v0  w0  ipx  ipy  <px  1 
Mid-side  Nodes  -  {w0  i px  yy) 

DISPLACEMENT  FIELD: 


.  u  =  u0  +  2^x  +  z2tx 

v  =  v0  +  z<Py  +  z24>y 

w  =  w0 

where ; 

u0,  v0,  <fx,  <f>y  :  {1  x  y  xy}T  {a} 

w0 >  v'x»  'Py  :  n  xy  x2  xy  y2  x3y  xy3}1  {3} 

STRESS  FIELD: 


i.  From  constitutive  relations  -  oj  =  C-jjc;j  (orthotropic  mat.) 

cxx  =  f(z2«  *2>  y2) 

Oyy  =  f(z2,  x2,  y2) 

oxy  =  f(z2,  x2,  y2) 

ii.  From  equilibrium  considerations  -  ojj.j  =  0 
oxz  =  f(z3,  x,  y) 

‘  OyZ  =  f(z3,  x,  y) 

ozz  =  f (z 3 ) 


QHD23 


NODAL  DEGREES  Or  FREEDOM 


{uo  v0  '•')  U'x  «y  $x  .yT 


DISPLACEMENT  FlfLi 


u  +  u0  +  + 


v  =  v0  +  z;'y  +  z2?y 


whore; 


V  V  V  -x*  V  V  4y,  : 


U  x  y  :<y}T{:,} 


STRESS  FIELD: 


i.  Frv...:  constitutive  relations  -  o-j  =  C -j j -  -j j  (orthotropic  mat.) 
‘-xx  s'T*(2r,  x,  y) 

:*yy  =  »  x’  y) 

cxy  ,  x,  y) 

i.  Fro-i  constitutive  considerations  -  e 3  0 

J  J 

°xz  =  f(z-) 

CyZ  f  (  Z  5  ) 

oZ2  -  constant 


APPENDIX  IB  -  MODIFIED  KIRCHHOFF  FORMULATION 

QD32 


NODAL  DEGREES  OF  FREEDOM: 


r  .  T 

Corner  Nodes  iw0  v0  w  y-  y-  yx  Yy) 
Mid-Side  Nodes  (w 

dn 


DISPLACEMENT  FIELD: 


w  =  f(x,  y) 
u  =  u0  -  z(j^  +  «xj 

v  ■  v0  -  z(^  ♦  ty) 


where : 


’  vo  *  x  *  y 


(1  *y  xy)Y:*} 

w  =  { 1  x  y  xz  xy  y2  x3  x2y  xy2  y3  x4  x3y  xy3  y4  x4y  xy4:T{3} 


STRESS  FIELD: 


i.  From  constitutive  relations  -  ij  (orthotropic  mat.) 


-'xx 


=  f(z,  x2,  y2) 


ayy  =  f(z,  x2,  y2) 
-  f(z,  X2,  y2) 


Jxy 


ii.  From  equilibrium  considerations  -  -  0 


oxz  =  f(z2,  x2,  y:‘ ) 
oyz  =  f(z2,  x2,  y2) 
°zz  =  f(z3,  X,  y) 


ill 


/ 


w  =  f(x,  y) 


U.u0-z(^*Yx) 

V  ’  V0  -  z  (  Ty  *  '•*) 


where; 


U  x.  y}T{a} 


u0’  v0’  X’  y 

w  :  U  x  y  x*  xy  y"  x3  x'"y  xy2  y3  x”  x3y  x2y2  xy3  y“}  •£: 


STRESS  FIELD 


i.  From  constitutive  relations  -  =  C-jjEy  (orthotropic  mat.) 

cxx  =  f(z,  x2  ,  y  ) 

cyy  =  f(z,  X2 ,  y2 ) 

Cxy  =  f(z,  X2  ,  T  ) 

ii.  From  equilibrium  considerations  -  ’ij»ij  ~  ® 
uxx  =  f(z:,  x.  y) 

OyZ  =  f(z2,  X,  y) 
c2Z  =  f(z3) 


/ 

/ 

/" 


1 


A 


APPENDIX  II 


HYBRID-STRESS  formulation 


DISPLACEMENT  FIELD: 


where; 


STRESS  FIELD: 


QLS32 

u  =  u0  +  z;x  +  z2<>x 
v  =  v0  +  Zfy  + 
w  =  wn  +  Z4>2 


^o»  vof  ^x»  ^y  *  z 9  *y  •  x  y  xy}  {a} 


°x  =  (Si  +S2  X  +  63y  +  34 xy)  +  z(85  +B6  x  +  g7y  +  88xy) 

+  z2  (e9  +  g10x  +  8ny  -  p64xy) 

ay  *  (S12  +  013  X  +  el4y  -  64xy)  +  z(6l5  +  s16  x  +  s17y  +  's,3  xy) 

+  Z2  ( 3 19  -£ltx  -  ivy  +  34xy) 

cxy  "l  '3  +  "  3  3n  + -j  3n  )x  +  (-S2  -  -5  810  -  ^  822  )y 

+  (-S9  -  8 19  )xy]  +  z  [S23  +  324  x  +  025y  +  (-1  b9-  1  ^  )*y] 

'  +  z2[5^  +  e2i  x  +  e^y  +  (p3g+  pg^xy] 
cxz  =  ^"h  "  t  "1^10  +  '322  )  +  34y  +  (-89^  819  )x] 

+  f(h2-  z2)[86  +  325  +  e6y  -  ^(89  +s 19 ) x ] 

+  ^(-h3-  z3)[,<10  -  p64y  +  622  +  p(e9+  6ig  )x] 

ayz  (-H  -  z)  [  -  8 21  +  8 11  )  +  (-89  -  8 19  )y  -  34x] 

+  2^'~  +  ^  +  3 17  +  S 18  x ] 

I  “l 

+  3  (-h3-  Z3)f62|  +  p(B9  +  3 19  )y  -  8n+  ~2s4x] 

*  (39  +  2,9)  [-(h  +  z)2-  (-?-h.3-  3h3z  +  Z'p  ( 3h^_4h3zj_z^)  i 


*»v\ 


APPENDIX  III  -  MASS  MATRIX  FORMULATION 

The  mass  matrix  for  elements  under  development  is  easily  arrived  at 
by  considering  kinetic  energy  in  the  form 

r 

T  =  \  1  p(u2  +  v2  +  w2)dV 

V 

where  u,  v  and  w  represent  displacements,  p  is  the  mass  density  and  the 
dot  superscript  denotes  velocity.  Defining  velocities  in  terms  of  element 
shape  functions  gives 

1  =  y  [L] 

V 

which  is  the  classical  form 
T  =  j  [A]  [M]\  A; 

The  element  mass  matrix  [M]  is,  therefore,  specified  as 

[M]  =  f  pff^.HNy]  +  (NV}[NV]  +  {NW}[NW]  dV 

J 

V 

Note  that  the  shape  functions  [N^]  involve  distance  from  the  mid-plane 
of  the  element  to  a  layer  denoted  by  Z  and,  therefore,  the  mass  matrix 
definition  provided  not  only  represents  mid-plane  inertial  effects  but 
also  rotary  inertia  as  well. 


J  c(Nu}[Nu]  t  {NV}[NV]  +  fNw}[Nw]  dV  {a} 


APPENOIX  IV 


LARGE  DISPLACEMENT  FORMULATION 


Based  on  Green's  Strain  Tensor,  the  following  procedure  is  utilized 
to  obtain  the  large  displacement  and  the  geometric  stiffness  matrices. 

Let  N  be  shape  functions  relating  displacements  at  any  point  in  the 
element  {5}  to  nodal  displacements  {A}  such  that 

{«}  •  [N]{a} 

Also  let  denote  those  shape  functions  associated  with  the  itfl 

displacement  field  (i  -  u,v,w)  and  ",j"  denotes  the  differentiation  with 
respect  to  the  jth  coordinate,  i.e.,  t—  where  x,  =  x,  x,  =  y  and  x,  =  z. 

,  i<jj .  ... 

Then,  the  strain  exx  given  by  J 


can  be  written  as 


£XX  =  ^u?x^  +  2  I^U’X^  .^U’X^  +  i^V»x'  ^V’x^  +  ^w»x)  {W|^> 


Similarly;  the  shear  strain  c  can  be  represented  by 

Cvy  -j{Nu,yiT  +  {Nv,x}T|{A}T  +  {A}T J{NU ,X}T{NU ,y } + {Nv,x)T{Nv,y} 


*  *  NW  ’ X *  [  (Ai 


The  strain  field  in  indicial  notation  is  expressed  by 

{cij}  =  |[[fNi,j}T  +  +  {A}T[{Nk,1.}T{Nk,j}]{A}} 


Then  the  incremental  representation  becomes 

{5eij}  =  ?[{{Ni»j}T  +  {NJ’i)T}^  'A}  +  {6A}T[{Nk>i}T{Nk,j}]{i} 
+  {A}T  [{Nk,i>T{Nk,j}]  {6A}] 

But  the  second  term  can  be  expressed  as 
{A}T[{Nk,j}T{Nk,i}]  {?.&) 


Thus  combining  terms 


^cij^  =  2  £ [  f|i  ’J ^  +  *Nj^i}TJ{£A}  McA}^[{Nk,i}^{Nk,j}  +{Nk,j}^fNk,i}]  { 6A>J 


Let 


[B0J  =  ♦  (Bj,()T] 

[BL]  =  ]-V:T  [;Nk,j)T!Nk,j)*(Nk,J)T  !Nk,|)]  - 

Then 


{MT[  Myy] 
{A}"^[MXy] 

{a)t[mzz] 

CA>T[Mxz] 

(A}T[Myz] 


( j }  -  [Bq]  1  +  f  BL]'. : ; 

where  [B0]  is  the  linear  component  and  [B|_]  is  the  large  displacement  component 
Having  the  definitions  for  [BQ]  and  [B^];  the  small  and  large  displacement 
matrices  [K0]  and  [KL]  are  represented  as 

W  =  /  CB0]T[D]CB0]dV 

V 

W  =  J  j  [BlJT[D][B0]  +  [BL]T[D][BL]  +  [B0]T[D][BL]  jdV 

V 


/\-4 

•y* 


The  geometric  stiffness  matrix  is  "’so  derived  from  [BL]  and  it  has  the 
followi.  ;  form 


■/( 


OxxCMxxl  +  °yyCMyy]  +  °zz ^MzzJ  +  °xy£Mxy]  +  °XZ[MXZ] 


Where  the  a's  are  the  stress  components  and  again  integration  is  on  a 
layer  by  layer  basis. 


* 

V 

m 


I 

r 


